%=========================================================================%
%-------------------Taylor-Galerkin 1D Burgers----------------------------%
%=========================================================================%

%=========================================================================%
% Simple implementation of the Taylor-Galerkin method for the 1D inviscid
% Burgers equation with transmissive boundary conditions. Uses a constant
% mesh sizing to demonstrate how the scheme translates into the well-known
% Lax-Wendroff formulation. This scheme requires CFL <=1, but existence of
% spurious oscillations near discontinuities and shocks implies CFL <= 0.8
% is a more realistic limit.
%=========================================================================%

%=========================================================================%
% The 1st step is to assume [u(:,n+1) = u(:,n)+dt*u,t+((dt^2)/2)*u,tt], in
% which u,t = du/dt. Now, from the conservation law, [u,t = -F,x], where
% F = F(u) is the flux function, [F = 0.5*u^2] in this case. The 2nd order
% term is then [u,tt = (AF,x),x], where A = dF/du, the Jacobian term. This 
% higher order term will lead to an added numerical dissipation with small
% value, that stabilizes the central method resulting from SG-FEM.
% Now, the true terms are approximated as [ug = Ni*ui] and [Fg = Ni*Fi]. In
% other words, take the sum of the product between shape functions and
% nodal values as the approximation, using both the variable u(x,t) and the
% flux F(u(x,t)). Requiring that the resulting residual is zero in the
% average weighted sense, and selecting w(x) = transpose(N(x)) as the
% weight leads to a standard Galerkin (SG-FEM) discretisation. The 2nd
% order term must of course be reduced to a weak form, where the resulting
% Neumann BCs are set to zero. Setting an equal mesh size and carrying out
% the full discretization procedure for a 2-element assembly leads to the
% original Lax-Wendroff scheme.
%=========================================================================%

clear;
clc;

% Domain and Time:

xl = 0; % Leftmost point
xr = 1; % Rightmost point
dx = 0.0005; % Mesh size
x = xl:dx:xr; % Equispaced Mesh
nx = length(x); % Grid points count
x0 = 0.5; % IC breakpoint
N = find(x == x0); % Breakpoint index

t0 = 0; % Initial time
tf = 4*0.3; % Final time
dt = 0.0003; % Timestep
t = t0:dt:tf; % Time points
nt = length(t); % Time point count

% Initial data and pre-alloc.:

u = zeros(nx,nt); % Data storage
u(1:N,1) = sin((pi*x(1:N))/x0); % Initial condition as Riemann Problem
u(N+1:nx,1) = 0;

% Solver:

k = 1; % Time counter start
while k ~= nt
    for i = 2:nx-1
        F = 0.5*u(:,k).^2; % Burgers Flux
        ND = -((dt^2)/(4*(dx^2)))*((-u(i-1,k)-u(i,k))*F(i-1)+(u(i-1,k)+...
            2*u(i,k)+u(i+1,k))*F(i)+(-u(i,k)-u(i+1,k))*F(i+1)); % Dissip.
        u(i,k+1) = u(i,k)-(dt/(2*dx))*(F(i+1)-F(i-1))+ND; % LW scheme
    end
    
    u(1,k+1) = u(2,k+1); % Left BC
    u(nx,k+1) = u(nx-1,k+1); % Right BC
    k = k+1; % Counter Adv.
    
end

% Post-processing:

figure(1)
hold on
subplot(1,2,1)
grid('minor')
plot(x,u(:,1:100:k),'k')
xlabel('x')
ylabel('u(x,t)')
title('Solution time-lapse')
subplot(1,2,2)
contourf(x,t,transpose(u),20),colormap('jet')
xlabel('x')
ylabel('t')
title('Characteristics plane')
hold off